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In this paper, we present a method to construct the eigenspace of the normal-state electrons 
moving in a 2D square lattice in presence of a perpendicular uniform magnetic field which imposes 
(quasi)-periodic boundary conditions for the wave functions in the magnetic unit cell. An exact 
unitary transformations are put forward to correlate the discrete eigenvectors of the 2D electrons 
with those of the Harper's equation. The cyclic-tridiagonal matrix associated with the Harper's 
equation is then tridiagonalized by another unitary transformation. The obtained eigenbasis is 
utilized to expand the Bogoliubov-de Gennes equations for the superconducting vortex lattice state, 
which showing the merit of our method in studying the large-sized system. To test our method, we 
have applied our results to study the vortex lattice state of an s-wave superconductor. 

PACS numbers: 

I. INTRODUCTION 

Vortex states of type-II superconductors has received greater attention in recent years. Theoretical formalism 
describing this effect is the Bogoliubov-de Gennes (BdG) approach^, which can be viewed as real-space extension of the 
Bardeen-Cooper-Schrieffer (BCS) theory. This method allows one to reveal effects of imperfections in superconductors, 
such as impurities, surfaces, as well as field-induced vortices which we are concerned with in this paper. In recent years 
there are numerous studies on the superconducting vortex lattice state by solving the discrete BdG equations and 
consequently diagonalization of the BdG mean-field Hamiltonian on a two-dimensional tight-binding lattice^ii^ii^i^i^. 
However, the size of the unit cell of the vortex lattice, which is inversely proportional to the amplitude of the magnetic 
field, is limited by computer resources since the dimension of the BdG equations grows with system size. Therefore 
early numerical works on small-size unit cells are limited to high magnetic fields over ten Tesla, which is stronger than 
used in most experiments, and no remarkable progress has been made over the past decade due to the time consuming 
of full diagonalization (i.e. all eigenvalues and eigenvectors) of the mean-field BdG Hamiltonian. In fact in BCS-type 
superconductors, electrons near the Fermi level bind into Cooper pairs by exchanging virtual bosons such as phonons, 
excitons or plasmons etc. Therefore, there exists an energy cutoff, which equals approximately the characteristic 
energy of the bosons such as the Debye phonon frequency of conventional superconductors, and correspondingly only 
the electronic states lying near the Fermi surface within an energy shell are necessary to be explored. For the vortex 
problem, the most appropriate starting point is to find the relevant electronic states which participate in BCS pairing 
and forming of the superconducting vortex lattice when an external magnetic field is applied. The eigenequation 
describing this states is a 2D difference equation formulated on a magnetic unit cell which is twice the size of that 
of the superconducting vortex lattice. This eigenvalue problem is also demanding when the system size is large even 
though only a truncated eigenspace is desired. 

In this paper we present an exact reduction of the Hermitian matrix associated with the 2D discrete equation into 
a tridiagonal matrix, which composed of two consecutive unitary transformations. First we reduces the 2D discrete 
equation that describes electrons moving in a magnetic unit cell into the famous Harper's equation. Algebraically, this 
unitary transformation reduces the Hermitian matrix into a cyclic tridiagonal matrix corresponding to the Harper's 
equation. Then by another exact transformation the cyclic tridiagonal matrix is further reduced into a tridiagonal 
form. The exact reduction greatly lessens the computational burden in numerical methods. Ultimately we diagonalize 
the tridiagonal matrix utilizing available standard software packages to find the appropriate eigenstates near the Fermi 
level, i.e. the truncated eigenspace, and then expand and diagonalize the BdG equations in this truncated eigenbasis. 

This paper is organized as follows. In Section |TT1 we derive the BdG equations expanded in terms of the truncated 
eigenbasis of the normal-state electrons in the magnetic field. The Hermitian matrix associated with the 2D tight- 
binding electrons on a 2D square lattice in a magnetic field is reduced into a tridiagonal form in Section IIIII In 
Section HVl the vortex lattice state of an s-wave superconductor is studied as a test of our method. Section fVl gives 
the concluding remarks. 
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II. THE BDG EQUATIONS FOR VORTEX LATTICE STATES 

In this work, we adopt a BCS-type mean-field Hamiltonian defined on a two-dimensional (2D) square lattice, 



ij.o- i,J 



where Ay = ^(ciiCjj^ — Ci|Cj|) for spin-singlet pairing''. In an external uniform magnetic field applied in the z-direction, 
the hopping integral acquires the Peierls phase factor as 



/.27r r \ j i^{m^,my), j ^ {m^ + l,my) 



Here t denotes the nearest-neighbor hopping integral. We choose the Landau gause with A — B{0,x,0) and the 
screening field induced by the supercurrent is neglected for extreme type-II superconductors. (j)a = h/eis the electronic 
flux quantum. Hereafter we use pair of integers i = {mx,my) as index of the site in the square lattice to denote the 
X and y coordinates. In the Nambu representation, the above Hamiltonian can be written as 



H 



h- fii A 



= *M A. r« f * (4) 
\ A* -h* + ^il J ^ ' 

where is the Nambu creation (annihilation) operator defined as ^'^ = {c\p ci|, • • • , c||, Ci|, • • • , c^|, cki) with 

K the total number of lattice sites, h and A are K x K matrices with elements {h)ij — and (A)ij — Ay, 
respectively. / is the K x K identity matrix. The mean-field Hamiltonian can be diagonalized by solving the following 
BdG equations, 

iy-^(5y Ay ^^ f U„(j) \ _ f U„(i) 



a;, -t;. + ^<5y J \ «„(j) J~^A vl{i) 1 



which can be viewed as Schrodinger-like equations for the electron and hole amplitudes of a BdG quasiparticle. The 
pairing potential Ay couples the u and v components and satisfies the self-consistent condition 

Ay=F ^"(iK(j)tanhf^), (6) 



|-B„|<Bd 



where Ed is the Debye-type cutoff energy of the pairing interaction. The BdG equations ([5]) can be expressed 
compactly in a matrix form 

h-fii A \ f u\ uy 



A* -h* + fil J \ v J yv 

with u and v i^-dimensional vectors. 

Abrikosov vortices, each of which carries one superconducting flux quantum $o = h/2e, are created and form a 
lattice structure in a type-II superconductor if one applies a magnetic field {Bd < B < Bc2)- The vortex lattice 
causes periodic modulation of the pairing potential and accordingly yields energy bands of BdG quaiparticles. To 
study this effect in our study we adopt the concept of magnetic unit cell (MUG) whose size is twice that of the unit 
cell of the vortex lattice and accordingly each MUG accommodates one electronic flux quantum c/jq = 2$o- Here for 
illumination of our method, we study the square vortex lattice which is aligned with the underlying crystalline lattice. 
The unit cell size of the vortex lattice is x , corresponding to a uniform magnetic field B = / (Nxo)'^ . Each 
MUG accommodates two adjacent vortices in the y direction. Therefore the MUG is of size x Ny with Ny = 2Nx- 
The whole system is composed of x My MUG's. Thus the whole system has M^MyN^Ny lattice sites. For later 
convenience, we introduce a dimensionless parameter a = Ba^/cjyQ — l/{NxNy) denoting the ratio of magnetic flux 
per plaquette to the electronic flux quantum 0o- 
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In the Abrikosov vortex lattice state, the BdG equations ([5]) is symmetric under magnetic translation with the 
translation vector R — l^N^e^ + lyNyBy. Due to this magnetic translational symmetry in the x and y direction, the 
quasiparticle amplitudes can be expressed in the magnetic Bloch form as 



^(i) 



(8) 



where the magnetic Bloch wave vector k = j^^j^^x ~r xf-^ 
reduces Eq. ([7|) to the new BdG equations for and 



■Gy with Ix^y = 0, 1, • • • , Mx^y — 1. This transformation 



h^- III 



-{k 



where the matrix elements of the k-dependent matrices and A"^ are {h^)ij — tijt 



-ik-(i-j) 



The quasiparticle amplitudes and V" 
direction 



(A'')ij 



satisfy the quasi-periodic boundary conditions with period N.^ 



Aije 



u^{mx + Nx,my) 
v^irrix + Nx,my) 



J2TrmyN^a k 



V {mx,my) 



(9) 

ik-(i-j) 

7 

along the x 
(10) 

(11) 



while they are periodic in the y direction with period Ny. The {mxTmyYs in Eqs. (|9ll0llip are restricted to sites 
within one MUG with rrix^y = 0, 1, • • • , N^^y — 1. The above procedure reduces the Hermitian matrix with linear 
dimension 2MxMyNxNy [ Eq. ([7]) ] into direct sum of M^My block matrices, each of which is labeled by k and has 
linear dimension 2NxNy [ Eq. ([9]) ]. For each quasimomentum k, Eq. ([9]) is diagonalized along with the boundary 
conditions and then the whole solutions of all k are used by the following equation 



J2 u^{i)[v^m*e'''-^'-'^tanh( 

\EY^\<Eo 



\2kBT 



(12) 



to achieve self-consistence. 

In the literature typical size of the unit cell of the vortex lattice studied by previous works was limited around 
20 X 202i2ii^i^i^. Such a small unit cell size corresponds to a magnetic field as large as i? = $o/(20a)^ w 32 Tesla, which 
is much higher than used by most experiments, if one assumes a typical lattice constant a w4A. Therefore one should 
find the way to diagonalize the BdG Hamiltonian (Hermitian matrix) with larger scale in order to match numerical 
calculation with experimental data. Although one can take advantage of the sparse nature of the BdG Hamiltonian 
Cl, we think that iterative methods, such as the Lanczos algorithm, are not appropriate for this problem because they 
are designed to compute a few eigenvalues (eigenvectors) with largest/smallest magnitudes. 

To study the vortex lattice state with larger unit cell and correspondingly weaker and realistic magnetic field, we re- 
express the real-space BdG equations Eq. © in the diagonal representation of which describes the 2D tight-binding 
electrons in presence of a magnetic field, 



^^^^ 



(13) 



where (f^ obeys the same boundary condition as Eq.dTU]). According to the BCS theory, only a fraction of electrons 
in the energy shell Ejj around the Fermi energy participate in the Cooper pairing. Therefore we should first get the 
eigenstates ip^ from Eq. (fT5|) with energies |e^ — /i| < E^ relative to the Fermi level, which will be addressed in Section 



mil Here the quasiparticle amplitudes and are expanded in the basis functions (p^ and {ifiq^)*, respectively. 



This reduces Eq. ^ to. 



E 



(4 



A'' 



anil) 



E'^ 



(14) 



(15) 



where the matrix element Ap ^ is calculated according to, 

A^^^ ^ (p^yA^i^-^y = 5][<^^(i)]*Aije-*-(--')[^-'^(j)]*, 



(16) 
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while from Eqs. (fT2|) and ((T4)) . we have 

Aij = ^ E ^p(i)^^^j)««(p)[&n(9rtanh('^') (17) 

The Eqs. (fT51) - P?)l are solved iteratively until self consistence is satisfied. Eventually we can calculated the local 
density of states, which is proportional to the differential tunneling conductance, from the energy spectrum and wave 
functions, 

p(i, E)^J2 - E^) + Iv^M'SiE + El). (18) 

k,n 

At the present stage, we have expressed the BdG equations in the truncated eigenbasis of . The issue now is how to 
compute this truncated eigenbasis, i.e. the eigenstates of h)^ lying within an energy shell Ed around the Fermi level. 

Utilizing standard computational algorithm^, it will be rather time-consuming to compute some selected eigenstates 
of a large matrix as h)^, whose size N x N grows fastly with the length scale of the MUG, by tridiagonalizing the 
matrix numerically. Even after taking advantage of the sparse nature of h^, we find that the iterative methods, 
such as the Lanczos algorithm, are not quite appropriate for this problem because they are most efficient for finding 
largest/smallest eigenvalues (eigenvectors). In the following sections, we solve this issue by showing that can be 
tridiagonalized exactly by two unitary transformations. Then we appeal to standard packages such as LAPACK^ to 
compute the desired eigenstates of the resulting tridiagonal matrix within an energy range. 



III. 2D TIGHT-BINDING ELECTRONS IN A MAGNETIC UNIT CELL 



In this section, we show in detail the method of reduce the matrix to a tridiagonal matrix exactly. The 
eigenequation of (p^ [Eq. (|13p ]. i.e. the discrete Schroodinger equation describing a 2D free electron moving in a 
perpendicular uniform magnetic field in a square lattice, can be written in an explicit form 

e'''^^^{m, + l,my)+e-'''-^^im,~l,my) + e'^^^"'^"+''y'>^^{m,,my + l) + 

(19) 

where = £n/(^0 ^^'^ obeys the quasi-periodic boundary condition along the x direction and periodic boundary 
condition along the y direction 



ip^{mx,my + Ny) = Lp^(mx,my). 



(20) 



First wc find that the eigenfunction ip^ is related to the cigcnfunction of the Harper's equation by a unitary 
transformation. Explicitly, 



Ar„-1 



1 - 



gl{m^ + lN^). 



(21) 



Substituting the above equation into Eq. p9p . one readily find that satisfies the Harper equation 



e''=^3j:(m+l) + e 



(m - 1) + 2cos(27rma + ky)g^^{m) = e^g^{m), 



(22) 



Here m = 0,1, ■ ■ ■ , N ~ I with N = N^Ny. g satisfies the periodic boundary condition gl{m + N) — gl{m). In the 
matrix form, the Harper equation can be expressed as 



(23) 



where 



(24) 



ajv-1 / 
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with ttm = 2 cos(27rma + fcy). And the eigenvector gjj = {g^{0), ■ ■ ■ , g^{N—l))'^. The Harper's equation can be viewed 
as discrete analogue of the Schrodinger equation of the one-dimensional quantum harmonic oscillator. Therefore, the 
energy eigenfunctions are discrete analogues of the Hermit-Gaussian-type wave functions. 

The periodic tridiagonal matrix can be further reduced to a tridiagonal matrix by another unitary transfor- 
mation. For simplicity, we only show the procedure for the ky = case and the following discussion can be readily 
generalized for ky ^ 0. The transformation of wave vectors from 17 to / is as follows. 



/(O) = .9(0), 
j^l^ ^ e''°"g(i)-l-' 

/(2) = 



'g(2)+e 



V2„ 



'g{N-i) 



/(f -1) = 
/(f) =. 9(f), 



/( 



V2 



= g(JV-2) 



I) 



V2 



(25) 



V2 



f{N-2) = 
f{N-l) = 



V2 



%/2 



Substituting the above relations into Eq. ([22|l . we have the eigenequation for /, which reads, 

Tnk /'k ~k /'k 

J n ^nJ n 

where T is an iV x tridiagonal matrix, 

/ flo V2 

ai 1 
1 • • 



(26) 



V2cos(f fc,) 



■\/2cos(-jfcj. 



i-\/2si 



svai^kx 



-iV2sin(f fc,) 

a jv 1 



1 



a N 



\ 



(27) 



1 ai / 



and/k = (/k(0),... ,/k(7V-l)). _ ....... 

After two consecutive unitary transformations, we have successfully reduce the Hermitian matrix h)^ into a tridi- 
agonal matrix T^. Here we emphasize that the reduction is exact without any numerical assumption and takes no 
CPU time compared with the numerical reduction. Then the eigenproblem of the tridiagonal matrix can be solved 
using standard packages such as LAPACK. 



IV. AN EXAMPLE: VORTEX LATTICE STATES OF A TYPE-II s-WAVE SUPERCONDUCTOR 



In this section, we illustrate how our method is applied in solving the BdG equation for the vortex lattice states 
of an s-wave superconductor. The microscopic parameters used in this paper are as follows. As a model calculation 
we set the relevant parameters as follows, /i — — 3t which gives rise to an almost circular Fermi surface with the 
Fermi wave vector k-p « 1.03a^^ and Fermi velocity v-p « l.^lta/fi. The on-site attractive interaction V = 2t. The 
Debye-type energy cutoff Eo = 0.5i. This set of parameters results in an s-wave pairing potential Aq ~ 0.065t in the 
zero-temperature limit with the estimated coherence length = Tiv-p / it « 9a. 

The model calculation is carried out for a system composed of x My = 40 x 20 MUC's with each MUG of size 
Nx X Ny = 80 X 160 which corresponds to a magnetic field B — (f)^/ {NxNyaY ~ 2.0 Tesla if the lattice constant is 
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set as 4A. Therefore, for each k of the total 800 quasimomenta, we employ standard LAPACK routine to diagonalize 
the 12800 X 12800 tridiagonal matrix fEq. [77]) and find that there are approximately 1173 eigenstates {/^} whose 
eigenenergies lying within the energy range je'' — A^l < Ed- We can obtain the eigenstates of the Harper's equation {g^} 
by the inverse transformation of Eq. (P5|) . Then substituting into Eq. (pi]) we successfully obtain the truncated 
eigenbasis {t/?5^}, in which the BdG equations (fT5|) are expressed as the 2 x 1173-dimensional eigenvalue problem 
and the matrix elements A^^ is calculated from Eq. (I16|) . After the BdG equations are diagonalized for each k, we 
substitute the quasiparticle amplitudes and b^^ into the self-consistent condition Eq. (|17p to compute the renewed 
values of the pairing potential. Eqs. p5|) - P?|) arc solved iteratively until convergence is reached. 
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FIG. 1: Spatial distribution of the magnitude of the s-wave pairing potential |A| in one magnetic unit cell of size 80 x 160. 
The X- and t/-axis are in unit of the lattice constant a. \A\ is in unit of the hoping integral t. 

In Fig. [1] we show the spatial variation of the self-consistent pair potential within one 80 x 160-sized magnetic 
unit cell, in which two superconducting vortices are situated. The s-wave pairing potential vanishes at the center of 
each of the two 80 x 80 squares and increases with the distance from the core center and recovers to its bulk value 
approximately with a length scale ^o- The variation of the pairing potential around the vortex core exhibits almost 
circular symmetry as shown in the figure. The reasons are twofold. Firstly, the Fermi level is far away from the van 
Hove singularity and accordingly the Fermi surface is approximately circular. Secondly, the impact from neighboring 
vortices which arranged squarely is weak because the distance between the adjacent vortices is about one order of 
magnitude larger compared to the the characteristic coherence length. 
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FIG. 2: Quasiparticle spectrum in the magnetic Brillouin zone. See text for details. 

Figure [2] displays the quasiparticle spectra along three high-symmetry lines in the magnetic Brillouin zone, where 
TX, XM and MF connect two of the three points: F (0,0), X = (#^,0) and M = {-j^, j^). As shown in 

the figure, the vortex bound states, which is localized in a isolated vortex line as revealed in RefsJ^iii, are broaden 
into energy bands in the superconducting vortex lattice owing to the interference effect. However due to the localized 
nature of the vortex states, the overlapping of the quasiparticle wave functions belonging to difference vortices is 
weak especially for the low-lying states. Consequently the bands with lower energies are flatter and the level spacing 
between pairs of the first few lowest-lying bands is of the order of Aq/Ep. These results are consistent with previous 
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worksii^ii^. 
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FIG. 3; LDOS as a function of energy at the vortex core center (solid line) and the inter-vortex site (dash-dotted line). 

In Fig. [3] we plot the local density of states (LDOS) as a function of energy at the vortex core center and the 
inter-vortex site. At the center of the vortex core, the LDOS is greatly enhanced at the energy approximately equal 
to Aq/Ef due to the strongly localized vortex bound states, while depressed around E — ±Ao as compared with the 
LDOS at the inter-vortex site. The model calculation shows the feasibility of our methods in studying the vortex 
lattice with large unit cell. 

V. CONCLUSION 

The discrete BdG equations developed in the 2D tight-binding lattice have been used to study the magnetic-field- 
induced superconducting vortex lattice state in the literature. The size of the system studied in previous works was 
limited due to the full diagonalization of the BdG hamiltonian directly. In this paper, we have extended this method 
by constructing a truncated eigenspace for the normal-state electrons moving on a 2D square lattice in presence of 
a uniform magnetic field. The motion of the electrons is governed by the vector potential, which impose a (quasi)- 
periodic boundary condition along the x and y directions of the magnetic unit cell. We have presented two consecutive 
unitary transformations to reduce the Hermitian matrix for the 2D electrons into a tridiagonal matrix exactly. By 
doing so, we have successfully related the desired eignbasis with that of the celebrated Harper's equation which is 
the eigenequation for a periodic-tridiagonal matrix. Then the second transformation is applied to further reduce the 
periodic tridigoanl matrix to a tridigonal one. This greatly reduces the cost of CPU time and helps us to treat systems 
with much larger size. To test our method and elucidate it more specifically, we have applied our results to study the 
vortex lattice states of an s-wave superconductor. The extension of our method to more sophisticated band structure 
as well as to 2D triangular or honeycomb lattice will be performed in future works. 
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